---
title: "MTurk data"
author: "Jeppe Sabroe Thegen"
date: '2022-06-23'
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F)
options(scipen = 999)
```

```{r}
#Load packages
library(tidyverse)
library(readxl)
library(psy)
library(stargazer)
library(ggeffects)
library(margins)
library(knitr)
library(tidycensus)
library(haven)
library(lmtest)
library(sandwich)

```

```{r Making Functions}
#Create function to normalize indexes between 0 and 1
range01 <- function(x, ...){(x - min(x, ...)) / (max(x, ...) - min(x, ...))}

#Create function to reverse variables
reverso <- function(x) {
y <- (x * (-1)) + max(x, na.rm = T)
return(y)
}

#Create a function to get the mode
getmode <- function(v) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}

force_bind = function(df1, df2) {
    colnames(df2) = colnames(df1)
    bind_rows(df1, df2)
}
```

```{r set wd}
setwd("~/Library/Mobile Documents/com~apple~CloudDocs/Speciale/Data")
```


```{r}
#Import MTurk data
MTurk <- read_csv("Psychological Distance to Climate Change_June 23, 2022_01.26.csv")
```



```{r}
#Initial Descriptions
##How many were blocked?
summary(MTurk$block %>% as.factor)
#Removing those
MTurk <- subset(MTurk, !(block %in% c(1, 2)))

#Where are they from?
summary(MTurk$countryName %>% as.factor)
#Removing those not from the US
MTurk <- subset(MTurk, (countryName %in% "United States"))

#How many did I show the wrong wording? 120
summary(MTurk$Q16 %>% as.factor)
#I code them as missing as answers are not meaningful
MTurk$Q16[MTurk$Q16 == 1] <- NA
MTurk$Q16[MTurk$Q16 == 2] <- NA
MTurk$Q16[MTurk$Q16 == 3] <- NA

#What is the distribution of responses?
summary(MTurk$Stimulus %>% as.factor)
#For four respondents, the stimulus was not recorded due to an error in the pilot phase.


```

```{r where are they?}
Byer <- read_excel("Byer.xlsx")
coords <- str_split(Byer$Coordinates, ",")
coords <- unlist(coords)
Byer$long <- coords[seq(1,length(coords),2)]
Byer$lat <- coords[-seq(1,length(coords),2)]
Byer$long <- Byer$long %>% as.numeric
Byer$lat <- Byer$lat %>% as.numeric

Byer <- subset(Byer, !(City %in% c("Brooklyn", "Manhattan", "Queens")))

usa <- map_data("state")
ggmap<- ggplot() + geom_polygon(data = usa, aes(x=long, y = lat, group = group),color= "darkslategrey", alpha = 0.6, size = 0.2, fill = "darkslategrey") + coord_fixed(1.3) + geom_point(data = MTurk, aes(x=LocationLongitude %>% as.numeric, y = LocationLatitude %>% as.numeric), size = 1, alpha = 0.3) + xlim(-125,-67) + ylim(25, 50) + theme(legend.position="bottom", axis.ticks.x=element_blank(), axis.ticks.y=element_blank(), panel.background = element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank()) + xlab("") + ylab("") + geom_point(data=Byer, aes(x=lat, y=long), color ="#E41A1C", alpha=0.8, shape = 15) + geom_text_repel(data=Byer, aes(x=lat, y = long, label = City), nudge_y = -0.4, color = "white", fontface = "bold", bg.colour = "black", bg.r = .1)

tikz(file = "map.tex",
     standAlone=F,
     width = 8,
     height = 4.5,
     sanitize =T,
     )

print(ggmap)

dev.off()

```

```{r}
#Cleaning it up

#Characteristics pre treatment
MTurk$age <- (MTurk$YOB %>% as.numeric-2022)*(-1)
MTurk$Party[MTurk$Party == "Prefer not to say"] <- NA
MTurk$republican <- ifelse(MTurk$Party=="Republican Party", 1, 0)
MTurk$male <- ifelse(MTurk$Gender=="Male", 1, 0)

#Makes education an ordered factor
MTurk$eduOrdered <- factor(MTurk$Education, ordered = TRUE, 
                                levels = c("Less than high school",
                                           "High school graduate",
                                           "Some college",
                                           "2 year degree",
                                           "4 year degree",
                                           "Professional degree",
                                           "Doctorate"))

#Replies post treatment
##Concern. Higher values indicate higher concern
MTurk$Q16[c(1:120)] <- NA #This was due to a survey error
MTurk$concernOverall <- MTurk$Q16
MTurk$concernOverall[MTurk$concernOverall=="Extremely worried"] <- 4
MTurk$concernOverall[MTurk$concernOverall=="Moderately worried"] <- 3
MTurk$concernOverall[MTurk$concernOverall=="Slightly worried"] <- 2
MTurk$concernOverall[MTurk$concernOverall=="Not worried at all"] <- 1
#It needs to be normalized before combining with the other items
MTurk$concernOverall <- range01(MTurk$concernOverall %>% as.numeric, na.rm=T)

############

MTurk$concernSocial <- MTurk$Q17_1
MTurk$concernSocial[MTurk$concernSocial=="Strongly disagree"] <- 7
MTurk$concernSocial[MTurk$concernSocial=="Disagree"] <- 6
MTurk$concernSocial[MTurk$concernSocial=="Somewhat disagree"] <- 5
MTurk$concernSocial[MTurk$concernSocial=="Neither agree nor disagree"] <- 4
MTurk$concernSocial[MTurk$concernSocial=="Somewhat agree"] <- 3
MTurk$concernSocial[MTurk$concernSocial=="Agree"] <- 2
MTurk$concernSocial[MTurk$concernSocial=="Strongly agree"] <- 1
MTurk$concernSocial <- range01(MTurk$concernSocial %>% as.numeric, na.rm=T)

MTurk$concernDistance <- MTurk$Q17_2
MTurk$concernDistance[MTurk$concernDistance=="Strongly disagree"] <- 7
MTurk$concernDistance[MTurk$concernDistance=="Disagree"] <- 6
MTurk$concernDistance[MTurk$concernDistance=="Somewhat disagree"] <- 5
MTurk$concernDistance[MTurk$concernDistance=="Neither agree nor disagree"] <- 4
MTurk$concernDistance[MTurk$concernDistance=="Somewhat agree"] <- 3
MTurk$concernDistance[MTurk$concernDistance=="Agree"] <- 2
MTurk$concernDistance[MTurk$concernDistance=="Strongly agree"] <- 1
MTurk$concernDistance <- range01(MTurk$concernDistance %>% as.numeric, na.rm=T)

MTurk$concernHypothetical <- MTurk$Q17_3
MTurk$concernHypothetical[MTurk$concernHypothetical=="Strongly disagree"] <- 7
MTurk$concernHypothetical[MTurk$concernHypothetical=="Disagree"] <- 6
MTurk$concernHypothetical[MTurk$concernHypothetical=="Somewhat disagree"] <- 5
MTurk$concernHypothetical[MTurk$concernHypothetical=="Neither agree nor disagree"] <- 4
MTurk$concernHypothetical[MTurk$concernHypothetical=="Somewhat agree"] <- 3
MTurk$concernHypothetical[MTurk$concernHypothetical=="Agree"] <- 2
MTurk$concernHypothetical[MTurk$concernHypothetical=="Strongly agree"] <- 1
MTurk$concernHypothetical <- range01(MTurk$concernHypothetical %>% as.numeric, na.rm=T)

MTurk$concernTime <- MTurk$Q17_4
MTurk$concernTime[MTurk$concernTime=="Strongly disagree"] <- 7
MTurk$concernTime[MTurk$concernTime=="Disagree"] <- 6
MTurk$concernTime[MTurk$concernTime=="Somewhat disagree"] <- 5
MTurk$concernTime[MTurk$concernTime=="Neither agree nor disagree"] <- 4
MTurk$concernTime[MTurk$concernTime=="Somewhat agree"] <- 3
MTurk$concernTime[MTurk$concernTime=="Agree"] <- 2
MTurk$concernTime[MTurk$concernTime=="Strongly agree"] <- 1
MTurk$concernTime <- range01(MTurk$concernTime %>% as.numeric, na.rm=T)

############Efficacy
MTurk$efficacyOthers <- MTurk$Q19_1
MTurk$efficacyOthers[MTurk$efficacyOthers=="Strongly disagree"] <- 1
MTurk$efficacyOthers[MTurk$efficacyOthers=="Disagree"] <- 2
MTurk$efficacyOthers[MTurk$efficacyOthers=="Somewhat disagree"] <- 3
MTurk$efficacyOthers[MTurk$efficacyOthers=="Neither agree nor disagree"] <- 4
MTurk$efficacyOthers[MTurk$efficacyOthers=="Somewhat agree"] <- 5
MTurk$efficacyOthers[MTurk$efficacyOthers=="Agree"] <- 6
MTurk$efficacyOthers[MTurk$efficacyOthers=="Strongly agree"] <- 7
MTurk$efficacyOthers <- range01(MTurk$efficacyOthers %>% as.numeric, na.rm=T)


MTurk$efficacySelf <- MTurk$Q19_2
MTurk$efficacySelf[MTurk$efficacySelf=="Strongly disagree"] <- 1
MTurk$efficacySelf[MTurk$efficacySelf=="Disagree"] <- 2
MTurk$efficacySelf[MTurk$efficacySelf=="Somewhat disagree"] <- 3
MTurk$efficacySelf[MTurk$efficacySelf=="Neither agree nor disagree"] <- 4
MTurk$efficacySelf[MTurk$efficacySelf=="Somewhat agree"] <- 5
MTurk$efficacySelf[MTurk$efficacySelf=="Agree"] <- 6
MTurk$efficacySelf[MTurk$efficacySelf=="Strongly agree"] <- 7
MTurk$efficacySelf <- range01(MTurk$efficacySelf %>% as.numeric, na.rm=T)

MTurk$efficacyExternal <- MTurk$Q19_3
MTurk$efficacyExternal[MTurk$efficacyExternal=="Strongly disagree"] <- 1
MTurk$efficacyExternal[MTurk$efficacyExternal=="Disagree"] <- 2
MTurk$efficacyExternal[MTurk$efficacyExternal=="Somewhat disagree"] <- 3
MTurk$efficacyExternal[MTurk$efficacyExternal=="Neither agree nor disagree"] <- 4
MTurk$efficacyExternal[MTurk$efficacyExternal=="Somewhat agree"] <- 5
MTurk$efficacyExternal[MTurk$efficacyExternal=="Agree"] <- 6
MTurk$efficacyExternal[MTurk$efficacyExternal=="Strongly agree"] <- 7
MTurk$efficacyExternal <- range01(MTurk$efficacyExternal %>% as.numeric, na.rm=T)

############Policy support

MTurk$hightax <- MTurk$Q20_1
MTurk$hightax[MTurk$hightax=="Strongly disagree"] <- 1
MTurk$hightax[MTurk$hightax=="Disagree"] <- 2
MTurk$hightax[MTurk$hightax=="Somewhat disagree"] <- 3
MTurk$hightax[MTurk$hightax=="Neither agree nor disagree"] <- 4
MTurk$hightax[MTurk$hightax=="Somewhat agree"] <- 5
MTurk$hightax[MTurk$hightax=="Agree"] <- 6
MTurk$hightax[MTurk$hightax=="Strongly agree"] <- 7
MTurk$hightax <- range01(MTurk$hightax %>% as.numeric, na.rm=T)

MTurk$airtravel <- MTurk$Q20_2
MTurk$airtravel[MTurk$airtravel=="Strongly disagree"] <- 1
MTurk$airtravel[MTurk$airtravel=="Disagree"] <- 2
MTurk$airtravel[MTurk$airtravel=="Somewhat disagree"] <- 3
MTurk$airtravel[MTurk$airtravel=="Neither agree nor disagree"] <- 4
MTurk$airtravel[MTurk$airtravel=="Somewhat agree"] <- 5
MTurk$airtravel[MTurk$airtravel=="Agree"] <- 6
MTurk$airtravel[MTurk$airtravel=="Strongly agree"] <- 7
MTurk$airtravel <- range01(MTurk$airtravel %>% as.numeric, na.rm=T)

MTurk$electriccars <- MTurk$Q20_3
MTurk$electriccars[MTurk$electriccars=="Strongly disagree"] <- 1
MTurk$electriccars[MTurk$electriccars=="Disagree"] <- 2
MTurk$electriccars[MTurk$electriccars=="Somewhat disagree"] <- 3
MTurk$electriccars[MTurk$electriccars=="Neither agree nor disagree"] <- 4
MTurk$electriccars[MTurk$electriccars=="Somewhat agree"] <- 5
MTurk$electriccars[MTurk$electriccars=="Agree"] <- 6
MTurk$electriccars[MTurk$electriccars=="Strongly agree"] <- 7
MTurk$electriccars <- range01(MTurk$electriccars %>% as.numeric, na.rm=T)

############Intentions
MTurk$footprint <- MTurk$Q21_1
MTurk$footprint[MTurk$footprint=="Strongly disagree"] <- 1
MTurk$footprint[MTurk$footprint=="Disagree"] <- 2
MTurk$footprint[MTurk$footprint=="Somewhat disagree"] <- 3
MTurk$footprint[MTurk$footprint=="Neither agree nor disagree"] <- 4
MTurk$footprint[MTurk$footprint=="Somewhat agree"] <- 5
MTurk$footprint[MTurk$footprint=="Agree"] <- 6
MTurk$footprint[MTurk$footprint=="Strongly agree"] <- 7
MTurk$footprint <- range01(MTurk$footprint %>% as.numeric, na.rm=T)

MTurk$transport <- MTurk$Q21_2
MTurk$transport[MTurk$transport=="Strongly disagree"] <- 1
MTurk$transport[MTurk$transport=="Disagree"] <- 2
MTurk$transport[MTurk$transport=="Somewhat disagree"] <- 3
MTurk$transport[MTurk$transport=="Neither agree nor disagree"] <- 4
MTurk$transport[MTurk$transport=="Somewhat agree"] <- 5
MTurk$transport[MTurk$transport=="Agree"] <- 6
MTurk$transport[MTurk$transport=="Strongly agree"] <- 7
MTurk$transport <- range01(MTurk$transport %>% as.numeric, na.rm=T)

MTurk$vacationing <- MTurk$Q21_3
MTurk$vacationing[MTurk$vacationing=="Strongly disagree"] <- 1
MTurk$vacationing[MTurk$vacationing=="Disagree"] <- 2
MTurk$vacationing[MTurk$vacationing=="Somewhat disagree"] <- 3
MTurk$vacationing[MTurk$vacationing=="Neither agree nor disagree"] <- 4
MTurk$vacationing[MTurk$vacationing=="Somewhat agree"] <- 5
MTurk$vacationing[MTurk$vacationing=="Agree"] <- 6
MTurk$vacationing[MTurk$vacationing=="Strongly agree"] <- 7
MTurk$vacationing <- range01(MTurk$vacationing %>% as.numeric, na.rm=T)

###Manipulation checks
###Higher values are far away
MTurk$ManipulateCloseness <- MTurk$Q25_1
MTurk$ManipulateCloseness[MTurk$ManipulateCloseness=="Does not describe well"] <- 5
MTurk$ManipulateCloseness[MTurk$ManipulateCloseness=="Describes slightly well"] <- 4
MTurk$ManipulateCloseness[MTurk$ManipulateCloseness=="Describes moderately well"] <- 3
MTurk$ManipulateCloseness[MTurk$ManipulateCloseness=="Describes very well"] <- 2
MTurk$ManipulateCloseness[MTurk$ManipulateCloseness=="Describes extremely well"] <- 1
MTurk$ManipulateCloseness <- range01(MTurk$ManipulateCloseness %>% as.numeric, na.rm=T)

MTurk$ManipulateSocial <- MTurk$Q25_2
MTurk$ManipulateSocial[MTurk$ManipulateSocial=="Does not describe well"] <- 5
MTurk$ManipulateSocial[MTurk$ManipulateSocial=="Describes slightly well"] <- 4
MTurk$ManipulateSocial[MTurk$ManipulateSocial=="Describes moderately well"] <- 3
MTurk$ManipulateSocial[MTurk$ManipulateSocial=="Describes very well"] <- 2
MTurk$ManipulateSocial[MTurk$ManipulateSocial=="Describes extremely well"] <- 1
MTurk$ManipulateSocial <- range01(MTurk$ManipulateSocial %>% as.numeric, na.rm=T)

MTurk$ManipulateLikelihood <- MTurk$Q25_3
MTurk$ManipulateLikelihood[MTurk$ManipulateLikelihood=="Does not describe well"] <- 5
MTurk$ManipulateLikelihood[MTurk$ManipulateLikelihood=="Describes slightly well"] <- 4
MTurk$ManipulateLikelihood[MTurk$ManipulateLikelihood=="Describes moderately well"] <- 3
MTurk$ManipulateLikelihood[MTurk$ManipulateLikelihood=="Describes very well"] <- 2
MTurk$ManipulateLikelihood[MTurk$ManipulateLikelihood=="Describes extremely well"] <- 1
MTurk$ManipulateLikelihood <- range01(MTurk$ManipulateLikelihood %>% as.numeric, na.rm=T)

MTurk$ManipulateTimeframe <- MTurk$Q25_4
MTurk$ManipulateTimeframe[MTurk$ManipulateTimeframe=="Does not describe well"] <- 5
MTurk$ManipulateTimeframe[MTurk$ManipulateTimeframe=="Describes slightly well"] <- 4
MTurk$ManipulateTimeframe[MTurk$ManipulateTimeframe=="Describes moderately well"] <- 3
MTurk$ManipulateTimeframe[MTurk$ManipulateTimeframe=="Describes very well"] <- 2
MTurk$ManipulateTimeframe[MTurk$ManipulateTimeframe=="Describes extremely well"] <- 1
MTurk$ManipulateTimeframe <- range01(MTurk$ManipulateTimeframe %>% as.numeric, na.rm=T)

```

```{r Doing factor analyses}

factanal(na.omit(MTurk[, c("concernOverall", "concernSocial", "concernDistance", "concernHypothetical", "concernTime")]), factors = 1)
#Just like in the Danish case, it seems that time is difficult to phrase.
cronbach(na.omit(MTurk[, c("concernOverall", "concernSocial", "concernDistance", "concernHypothetical", "concernTime")]))

factanal(na.omit(MTurk[, c("efficacyOthers", "efficacySelf", "efficacyExternal")]), factors = 1)
cronbach(na.omit(MTurk[, c("efficacyOthers", "efficacySelf", "efficacyExternal")]))

factanal(na.omit(MTurk[, c("hightax", "airtravel", "electriccars")]), factors = 1)
cronbach(na.omit(MTurk[, c("hightax", "airtravel", "electriccars")]))

factanal(na.omit(MTurk[, c("footprint", "transport", "vacationing")]), factors = 1)
cronbach(na.omit(MTurk[, c("footprint", "transport", "vacationing")]))



```


```{r Creating indexes}

MTurk <- mutate(MTurk, concern_index = rowMeans(dplyr::select(MTurk, c("concernOverall","concernSocial","concernDistance", "concernHypothetical","concernTime")), na.rm = T))

MTurk <- mutate(MTurk, efficacy_index = rowMeans(dplyr::select(MTurk, c("efficacyOthers","efficacySelf", "efficacyExternal")), na.rm = T))

MTurk <- mutate(MTurk, polsupport_index = rowMeans(dplyr::select(MTurk, c("hightax","airtravel", "electriccars")), na.rm = T))
#Interestingly, this shows some difference from what I expected.

MTurk <- mutate(MTurk, intentions_index = rowMeans(dplyr::select(MTurk, c("footprint","transport", "vacationing")), na.rm = T))

```

```{r Showing distribution of indexes}

gg_concern_index <- ggplot(MTurk, aes(x=concern_index)) + geom_density() + geom_histogram(aes(y=..density..), alpha=0.5, bins = length(unique(MTurk$concern_index))/2+1) + geom_vline(aes(xintercept=mean(concern_index, na.rm=T)), linetype = "dashed") + theme_minimal()+ labs(y="Tæthed", x="Bekymring")

gg_efficacy_index <- ggplot(MTurk, aes(x=efficacy_index)) + geom_density() + geom_histogram(aes(y=..density..), alpha=0.5, bins = length(unique(MTurk$efficacy_index))-4) + geom_vline(aes(xintercept=mean(efficacy_index, na.rm=T)), linetype = "dashed") + theme_minimal()+ labs(y="Tæthed", x="Efficacy")

gg_polsupport_index <- ggplot(MTurk, aes(x=polsupport_index)) + geom_density() + geom_histogram(aes(y=..density..), alpha=0.5, bins=length(unique(MTurk$efficacy_index))-4) + geom_vline(aes(xintercept=mean(polsupport_index, na.rm=T)), linetype = "dashed") + theme_minimal()+ labs(y="Tæthed", x="Policy support")

gg_intentions_index <- ggplot(MTurk, aes(x=intentions_index)) + geom_density() + geom_histogram(aes(y=..density..), alpha=0.5, bins=length(unique(MTurk$intentions_index))-4) + geom_vline(aes(xintercept=mean(intentions_index, na.rm=T)), linetype = "dashed") + theme_minimal() + labs(y="Tæthed", x="Intentioner")

tikz(file = "TurkIndekser.tex",
     standAlone=F,
     width = 6,
     height = 4,
     sanitize =T,
     )

print(pg_MTurk_indekser <- plot_grid(gg_concern_index, gg_efficacy_index, gg_polsupport_index, gg_intentions_index, 
          ncol =2, align = "hv"))

dev.off()

```


```{r Analyses}

#####Concern
summary(lm(concern_index ~ Stimulus + Stimulus*eduOrdered %>% as.numeric, data = MTurk))
#Interessant^
summary(lm(efficacy_index ~ Stimulus + Stimulus*eduOrdered %>% as.numeric, data = MTurk))
summary(lm(polsupport_index ~ Stimulus + Stimulus*eduOrdered %>% as.numeric, data = MTurk))
summary(lm(intentions_index ~ Stimulus + Stimulus*eduOrdered %>% as.numeric, data = MTurk))



concern_m1 <- lm(concern_index ~ Stimulus, data = MTurk)
HC_concern_m1 <- coeftest(concern_m1,vcov=vcovHC(concern_m1))


concern_m2 <- lm(concern_index ~ Stimulus + male + age + republican+ eduOrdered , data = MTurk)
HC_concern_m2 <- coeftest(concern_m2,vcov=vcovHC(concern_m2))

stargazer(concern_m1, concern_m2, type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne, MTurk",
          header = F,
          single.row = T,
          se = list(HC_concern_m1[,"Std. Error"],HC_concern_m2[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = append(c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid", "Mand", "Alder", "Republikaner"),levels(MTurk$eduOrdered)[-1]),
          out = "MTurkbekymring.tex",
          label="MTurkbekymring",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

concern_m3 <- lm(concern_index ~ Stimulus + male, data = MTurk)
HC_concern_m3 <- coeftest(concern_m3,vcov=vcovHC(concern_m3))

concern_m4 <- lm(concern_index ~ Stimulus  + Stimulus*male, data = MTurk)
HC_concern_m4 <- coeftest(concern_m4,vcov=vcovHC(concern_m4))


stargazer(concern_m1, concern_m3, concern_m4, type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne (betinget af køn)",
          header = F,
          single.row = T,
          se = list(HC_concern_m1[,"Std. Error"],HC_concern_m3[,"Std. Error"], HC_concern_m4[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid", "Mand", "Geografisk × Mand", "Hypotetisk × Mand", "Social × Mand", "Tid × Mand"),
          out = "MTurkbekymringbetingetafkon.tex",
          label="MTurkbekymringbetingetafkon",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
#Looks like there is a positive effect of time here (for women)

concern_m1_repub <- lm(concern_index ~ Stimulus + republican, data = MTurk)
HC_concern_m1_repub <- coeftest(concern_m1_repub,vcov=vcovHC(concern_m1_repub))

concern_m2_repub <- lm(concern_index ~ Stimulus + republican*Stimulus, data = MTurk)
HC_concern_m2_repub <- coeftest(concern_m2_repub,vcov=vcovHC(concern_m2_repub))

stargazer(concern_m1, concern_m1_repub, concern_m2_repub, type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne (betinget af parti)",
          header = F,
          single.row = T,
          se = list(HC_concern_m1[,"Std. Error"],HC_concern_m1_repub[,"Std. Error"], HC_concern_m2_repub[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid", "Republikaner", "Geografisk × Republikaner", "Hypotetisk × Republikaner", "Social × Republikaner", "Tid × Republikaner"),
          out = "MTurkbekymringbetingetafparti.tex",
          label="Mturkbekymringbetingetafparti",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

#####AGE####

concern_m1_age <- lm(concern_index ~ Stimulus + age, data = MTurk)
HC_concern_m1_age <- coeftest(concern_m1_age,vcov=vcovHC(concern_m1_age))

concern_m2_age <- lm(concern_index ~ Stimulus + age*Stimulus, data = MTurk)
HC_concern_m2_age <- coeftest(concern_m2_age,vcov=vcovHC(concern_m2_age))

stargazer(concern_m1, concern_m1_age, concern_m2_age, type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne (betinget af alder)",
          header = F,
          single.row = T,
          se = list(HC_concern_m1[,"Std. Error"],HC_concern_m1_age[,"Std. Error"], HC_concern_m2_age[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid", "Alder", "Geografisk × Alder", "Hypotetisk × Alder", "Social × Alder", "Tid × Alder"),
          out = "MTurkbekymringbetingetafalder.tex",
          label="Mturkbekymringbetingetafalder",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

concern_m1_edu <- lm(concern_index ~ Stimulus + eduOrdered, data = MTurk)
HC_concern_m1_edu <- coeftest(concern_m1_edu,vcov=vcovHC(concern_m1_edu))



stargazer(concern_m1, concern_m1_edu, type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på bekymring for klimaforandringerne (betinget af uddannelse)",
          header = F,
          single.row = T,
          se = list(HC_concern_m1[,"Std. Error"],HC_concern_m1_edu[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = append(c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid"),levels(MTurk$eduOrdered)[-1]),
          out = "MTurkbekymringbetingetafuddannelse.tex",
          label="Mturkbekymringbetingetafuddannelse",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

########

mConcernSocial <- lm(concernSocial ~ Stimulus, data = MTurk)
HC_mConcernSocial <- coeftest(mConcernSocial,vcov=vcovHC(mConcernSocial))
pvals_mConcernSocial <- list(HC_mConcernSocial[, 4])

mConcernGeo <- lm(concernDistance ~ Stimulus, data = MTurk)
HC_mConcernGeo <- coeftest(mConcernGeo,vcov=vcovHC(mConcernGeo))
pvals_mConcernGeo <- list(HC_mConcernGeo[, 4])

mConcernTime <- lm(concernTime ~ Stimulus, data = MTurk)
HC_mConcernTime <- coeftest(mConcernTime,vcov=vcovHC(mConcernTime))
pvals_mConcernTime <- list(HC_mConcernTime[, 4])

mConcernHypo <- lm(concernHypothetical ~ Stimulus, data = MTurk)
HC_mConcernHypo <- coeftest(mConcernHypo,vcov=vcovHC(mConcernHypo))
pvals_mConcernHypo <- list(HC_mConcernHypo[, 4])

stargazer(mConcernGeo, mConcernHypo, mConcernSocial, mConcernTime, 
          type="latex", out = "MTurk_underdimensioner_bekymring.tex",
          se = list(HC_mConcernGeo[,"Std. Error"],HC_mConcernHypo[,"Std. Error"], HC_mConcernSocial[,"Std. Error"],HC_mConcernTime[,"Std. Error"]),
          p=list(HC_mConcernGeo[,"Pr(>|t|)"],HC_mConcernHypo[,"Pr(>|t|)"],HC_mConcernSocial[,"Pr(>|t|)"],HC_mConcernTime[,"Pr(>|t|)"]),
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          header = F,
          single.row = T,
          intercept.bottom = FALSE,
          covariate.labels = c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid"),
          dep.var.labels = c("Geografisk", "Hypotetisk", "Social", "Tid"),
          dep.var.caption  = "Bekymringsmæssig underdimension",
          title = "Effekt af psykologiske afstande på bekymringsmæssige underdimensioner",
          label="bekymringerMTurk",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

#efficacy_m1 <- lm(efficacy_index ~ Stimulus + Stimulus*republican, data = MTurk[MTurk$republican %in% "0",])
efficacy_m1 <- lm(efficacy_index ~ Stimulus, data = MTurk)
HC_efficacy_m1 <- coeftest(efficacy_m1,vcov=vcovHC(efficacy_m1))

efficacy_m2 <- lm(efficacy_index ~ Stimulus + republican + age + male + eduOrdered, data = MTurk)
efficacy_m3 <- lm(efficacy_index ~ Stimulus + Stimulus*republican + male + age + eduOrdered, data = MTurk)
HC_efficacy_m3 <- coeftest(efficacy_m3,vcov=vcovHC(efficacy_m3))
HC_efficacy_m2 <- coeftest(efficacy_m2,vcov=vcovHC(efficacy_m2))


summary(efficacy_m3)
#Borderline significant effect of hypothetical distance on efficacy (lowers)
#It seems that democrats get *significantly less efficacious when hypothetical distance is high. Republicans are already not that efficacious. 
#Output for final
stargazer(m1_efficacy_index, m2_efficacy_index,efficacy_m1, efficacy_m2,type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_m1_efficacy_index[,"Std. Error"],HC_m2_efficacy_index[,"Std. Error"],HC_efficacy_m1[,"Std. Error"],HC_efficacy_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på følelse af efficacy",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tidsmæssig", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]),
          out = "MTurkefficacyalder.tex",
          label="MTurkefficacyalder",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
          column.labels   = c("DR Panelet", "MTurk"),
          column.separate = c(2, 2))


#13 har ikke ønsket at tilkendegive deres partitilhørsforhold.

efficacy_margs <- margins(efficacy_m3, at = list(republican = c(0,1)))
dfmargs_efficacy <- data.frame(summary(efficacy_margs))

ggmarg_efficacy <- ggplot(dfmargs_efficacy[dfmargs_efficacy$factor %in% "StimulusHypothetical",], aes(x = republican %>% as.factor, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Ikke-republikanere", "Republikanere")) + labs(x = "", y = "Marginal effect of hypothetical distance on efficacy",caption = "95 % confidence intervals")+ theme(plot.caption = element_text(hjust = 0.5), panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank())

tikz(file = "ggmarg_efficacy.tex",
     standAlone=F,
     width = 4,
     height = 4,
     sanitize =T,
     )

print(ggmarg_efficacy)

dev.off()

#Output models
polsup_m1 <- lm(polsupport_index ~ Stimulus, data = MTurk)
polsup_m2 <- lm(polsupport_index ~ Stimulus + republican , data = MTurk)
polsup_m3 <- lm(polsupport_index ~ Stimulus + republican + Stimulus*republican, data = MTurk)
stargazer(polsup_m1, polsup_m2, polsup_m3, type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001"),
          notes.append=FALSE)
summary(polsup_m1)
summary(margins(polsup_m3))
#Borderline significant effect of geographical on policy support (increases)
#There is a negative interaction effect of being republican.

pred.republican=data.frame(Stimulus="Social", republican=1)
pred.notrepublican=data.frame(Stimulus="Social", republican=0)
pred.republican.control=data.frame(Stimulus="Control", republican=1)
pred.notrepublican.control=data.frame(Stimulus="Control", republican=0)

predict(polsup_m3, newdata=pred.republican)
predict(polsup_m3, newdata=pred.notrepublican)
predict(polsup_m3, newdata=pred.republican.control)
predict(polsup_m3, newdata=pred.notrepublican.control)

#This function is equivalent to dydx in Stata
polsup_margs <- margins(polsup_m3, at = list(republican = c(0,1)))
dfmargs1 <- data.frame(summary(polsup_margs))

ggmargs_polsup_geo <- ggplot(dfmargs1[dfmargs1$factor %in% "StimulusGeographical",], aes(x = republican %>% as.factor, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Ikke-republikanere", "Republikanere")) + labs(x = "", y = "Geografisk afstand",caption = "") + theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank()) + coord_cartesian(ylim = c(-.25,.25))
#You can also do it with interplot
#interplot(m = polsup_m1, var1 = "StimulusSocial", var2 = "republican")

ggmargs_polsup_social <- ggplot(dfmargs1[dfmargs1$factor %in% "StimulusSocial",], aes(x = republican %>% as.factor, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Ikke-republikanere", "Republikanere")) + labs(x = "", y = "Social afstand",caption = "95%'s konfidensintervaller") + theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank()) + coord_cartesian(ylim = c(-.25,.25))

tikz(file = "ggmarg_polsup.tex",
     standAlone=F,
     width = 6,
     height = 4,
     sanitize =T,
     )

print(plot_grid(ggmargs_polsup_geo, ggmargs_polsup_social, 
          ncol =2))

dev.off()

#Calculate margins for social distance where significant
#-0.06697487 + (1.02602 * 0.06527658) I solved this in Wolfram Alpha

#From a Z-table I can see that it is only significant at the 69% level.
#So how many would be needed to see that effect?
#Input in Wolfram Alpha: 1926.8 respondents at 95% level
#-0.06697487 + (1.96 * (1.49994 /(sqrt(x)))) = 0
#1.49994(SD) comes from SD = 0.06527658(SE) * sqrt(528)

ggplot(dfmargs1[dfmargs1$factor %in% "StimulusHypothetical",], aes(x = republican, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Ikke-republikanere", "Republikanere")) + labs(x = "", y = "Marginal effect of social distance on policy support",caption = "95 \\% confidence intervals") + theme(plot.caption = element_text(hjust = 0.5), panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank())


# effects <- ggeffect(polsup_m1, terms = c("Stimulus", "republican"), 
#   ci.lvl = 0.95)
# ggplot(mydf, aes(x = x, y = predicted, colour = group)) +
#   geom_line()

intentions_m1 <- lm(intentions_index ~ Stimulus, data = MTurk)
intentions_m2 <- lm(intentions_index ~ Stimulus + republican, data = MTurk)
intentions_m3 <- lm(intentions_index ~ Stimulus + Stimulus*republican, data = MTurk)

stargazer(intentions_m1, intentions_m2, intentions_m3, type="text",
          star.char = c("†", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001"),
          notes.append=FALSE)

#Very little happens

```


```{r Margsconcern}

concernmargs <- margins(concern_m4, at = list(male = c(0,1)))
dfmargconcern1 <- data.frame(summary(concernmargs))

ggmargs_concernmargs <- ggplot(dfmargconcern1[dfmargconcern1$factor %in% "StimulusTime",], aes(x = male %>% as.factor, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Kvinder", "Mænd")) + labs(x = "", y = "", caption = "95 %'s konfidensintervaller", title = "") + theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank()) + coord_cartesian(ylim = c(-.25,.25))

tikz(file = "timemargsconcern.tex",
     standAlone=F,
     width = 4,
     height = 4,
     sanitize =T,
     )

print(ggmargs_concernmargs)

dev.off()


```

```{r Model assumptions}

container$res <- resid(m1_concern_index1)
qqnorm(res)
qqline(res)

plot(density(res))

ggplot(drpanel, aes(sample=res))+stat_qq()

```

```{r Political support comparison}

MTurk_polsup_m1 <- lm(polsupport_index ~ Stimulus, data = MTurk)
MTurk_polsup_m2 <- lm(polsupport_index ~ Stimulus + age + male + republican + eduOrdered, data = MTurk)
MTurk_polsup_m3 <- lm(polsupport_index ~ Stimulus + age + male + republican*Stimulus+ eduOrdered, data = MTurk)
MTurk_polsup_m4 <- lm(polsupport_index ~ Stimulus + age + male*Stimulus + republican+ eduOrdered, data = MTurk)
MTurk_polsup_m4 <- lm(polsupport_index ~ Stimulus + ageStimulus + male + republican+ eduOrdered, data = MTurk)


DR_polsup_m1 <- lm(support_index ~ Stimulus, data = drpanel)
DR_polsup_m2 <- lm(support_index ~ Stimulus + age + male + republican +eduyears + income, data = drpanel)

HC_MTurk_polsup_m1 <- coeftest(MTurk_polsup_m1,vcov=vcovHC(MTurk_polsup_m1))
HC_MTurk_polsup_m2 <- coeftest(MTurk_polsup_m2,vcov=vcovHC(MTurk_polsup_m2))
HC_MTurk_polsup_m3 <- coeftest(MTurk_polsup_m3,vcov=vcovHC(MTurk_polsup_m3))
HC_MTurk_polsup_m4 <- coeftest(MTurk_polsup_m4,vcov=vcovHC(MTurk_polsup_m4))
HC_DR_polsup_m1 <- coeftest(DR_polsup_m1,vcov=vcovHC(DR_polsup_m1))
HC_DR_polsup_m2 <- coeftest(DR_polsup_m2,vcov=vcovHC(DR_polsup_m2))

stargazer(DR_polsup_m1, DR_polsup_m2,MTurk_polsup_m1, MTurk_polsup_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Referencegruppe for uddannelse er 'Less than High School", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_polsup_m1[,"Std. Error"],HC_DR_polsup_m2[,"Std. Error"],HC_MTurk_polsup_m1[,"Std. Error"],HC_MTurk_polsup_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på støtte til politisk handling",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Politisk støtte",
          dep.var.labels = "Politisk støtte",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tidsmæssig", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "comp_polsup.tex",
          label="comp_polsup",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
           column.labels   = c("DR Panelet", "MTurk"),
          column.separate = c(2, 2))

stargazer(MTurk_polsup_m1, MTurk_polsup_m2,MTurk_polsup_m3,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for uddannelse er 'Less than High School", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_MTurk_polsup_m1[,"Std. Error"],HC_MTurk_polsup_m2[,"Std. Error"],HC_MTurk_polsup_m3[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på støtte til politisk handling, MTurk",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Politisk støtte",
          dep.var.labels = "Politisk støtte",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Hypotetisk", "Social","Tidsmæssig")), c("Alder", "Køn, Mand", "Republikaner")), levels(MTurk$eduOrdered)[-1]),c("Geografisk × Republikaner","Hypotetisk × Republikaner","Social × Republikaner", "Tid × Republikaner")),
          out = "polsup_repub.tex",
          label="polsup_repub",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
          column.labels=c("(3)", "(4)", "(5)"),
          model.numbers = F)

stargazer(MTurk_polsup_m1, MTurk_polsup_m2,MTurk_polsup_m4,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for uddannelse er 'Less than High School", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_MTurk_polsup_m1[,"Std. Error"],HC_MTurk_polsup_m2[,"Std. Error"],HC_MTurk_polsup_m4[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på støtte til politisk handling betinget af køn, MTurk",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Politisk støtte",
          dep.var.labels = "Politisk støtte",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Hypotetisk", "Social","Tidsmæssig")), c("Alder", "Køn, Mand", "Republikaner")), levels(MTurk$eduOrdered)[-1]),c("Geografisk × Mand","Hypotetisk × Mand","Social × Mand", "Tid × Mand")),
          out = "polsup_male.tex",
          label="polsup_male",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
          column.labels=c("(3)", "(4)", "(6)"),
          model.numbers = F)

```



```{r Intentions comparison}

MTurk_ints_m1 <- lm(intentions_index ~ Stimulus, data = MTurk)
MTurk_ints_m2 <- lm(intentions_index ~ Stimulus + age + male + republican + eduOrdered, data = MTurk)
MTurk_ints_m3 <- lm(intentions_index ~ Stimulus + age + male + republican*Stimulus+ eduOrdered, data = MTurk)
MTurk_ints_m4 <- lm(intentions_index ~ Stimulus + age + male*Stimulus + republican+ eduOrdered, data = MTurk)
MTurk_ints_m4 <- lm(intentions_index ~ Stimulus + age*Stimulus + male + republican+ eduOrdered, data = MTurk)


DR_ints_m1 <- lm(intention_index ~ Stimulus, data = drpanel)
DR_ints_m2 <- lm(intention_index ~ Stimulus + age + male + republican +eduyears + income, data = drpanel)

HC_MTurk_ints_m1 <- coeftest(MTurk_ints_m1,vcov=vcovHC(MTurk_ints_m1))
HC_MTurk_ints_m2 <- coeftest(MTurk_ints_m2,vcov=vcovHC(MTurk_ints_m2))
HC_MTurk_ints_m3 <- coeftest(MTurk_ints_m3,vcov=vcovHC(MTurk_ints_m3))
HC_MTurk_ints_m4 <- coeftest(MTurk_ints_m4,vcov=vcovHC(MTurk_ints_m4))
HC_DR_ints_m1 <- coeftest(DR_ints_m1,vcov=vcovHC(DR_ints_m1))
HC_DR_ints_m2 <- coeftest(DR_ints_m2,vcov=vcovHC(DR_ints_m2))

stargazer(DR_ints_m1, DR_ints_m2,MTurk_ints_m1, MTurk_ints_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Referencegruppe for uddannelse er 'Less than High School", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_ints_m1[,"Std. Error"],HC_DR_ints_m2[,"Std. Error"],HC_MTurk_ints_m1[,"Std. Error"],HC_MTurk_ints_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på handlingsintentioner",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Handlingsintentioner",
          dep.var.labels = "Handlingsintentioner",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tidsmæssig", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "comp_intentions.tex",
          label="comp_intentions",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
           column.labels   = c("DR Panelet", "MTurk"),
          column.separate = c(2, 2))
```

```{r Scepticism no comparison}

DR_scep_m1 <- lm(scepticism_index ~ Stimulus, data = drpanel)
DR_scep_m2 <- lm(scepticism_index ~ Stimulus + age + male + republican +eduyears + income, data = drpanel)

HC_DR_scep_m1 <- coeftest(DR_scep_m1,vcov=vcovHC(DR_scep_m1))
HC_DR_scep_m2 <- coeftest(DR_scep_m2,vcov=vcovHC(DR_scep_m2))

stargazer(DR_scep_m1, DR_scep_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_ints_m1[,"Std. Error"],HC_DR_ints_m2[,"Std. Error"],HC_MTurk_ints_m1[,"Std. Error"],HC_MTurk_ints_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på scepticisme",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Scepticisme",
          dep.var.labels = "Scepticisme",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tidsmæssig", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "scepticism.tex",
          label="scepticism",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r Efficacy hope}

DR_effhope_m1 <- lm(efficacy_index ~ treatment_hope, data = drpanel)
DR_effhope_m2 <- lm(efficacy_index ~ treatment_hope + age + male + republican +eduyears + income, data = drpanel)
DR_effhope_m3 <- lm(efficacy_index ~ treatment_hope + age + male + republican +eduyears, data = drpanel)

HC_DR_effhope_m1 <- coeftest(DR_effhope_m1,vcov=vcovHC(DR_effhope_m1))
HC_DR_effhope_m2 <- coeftest(DR_effhope_m2,vcov=vcovHC(DR_effhope_m2))
HC_DR_effhope_m3 <- coeftest(DR_effhope_m3,vcov=vcovHC(DR_effhope_m3))


######Laver en graf midt i det hele######
hope_efficacy <- drpanel %>%
  group_by(treatment_hope) %>%
  summarize(AvgEff = mean(efficacy_index, na.rm = T),
            SDEff = sd(efficacy_index, na.rm = T)
            ) %>% cbind(drpanel %>% count(treatment_hope))

hope_efficacy$low <- hope_efficacy$AvgEff - (1.96 * (hope_efficacy$SDEff/(sqrt(hope_efficacy$n))))

hope_efficacy$high <- hope_efficacy$AvgEff + (1.96 * (hope_efficacy$SDEff/(sqrt(hope_efficacy$n))))

hope_efficacy <- hope_efficacy %>% dplyr::select(-4)

gghope_efficacy<- ggplot(hope_efficacy, aes(x = AvgEff, y = treatment_hope)) + geom_point() +
geom_errorbar(aes(xmin=low, xmax = high), alpha=0.5, linetype = 3) +
  coord_cartesian(xlim = c(.55, .75), clip = "off") +
  labs(x = "", y = "", color = "", title = "Gennemsnitlig efficacy per gruppe, DR Panelet", caption = "95% konfidensintervaller vist") +
  theme_minimal() +
  theme(legend.position = "top", panel.grid.minor.y = element_blank(),panel.grid.major.y = element_blank())  + 
 # theme(plot.margin = unit(c(1, 1, 3, 1), "lines")) +
  annotate(geom = "text", x = c(.55, .75), y = -Inf, label = c("Lav efficacy", "Høj efficacy"), size = 2, color = "Dark Grey") + guides(fill = "none") + scale_y_discrete(labels = c("Ingen Treatment" = "Ingen stimulus","Hope" = "Håb", "None" = "Intet håb"))

######

stargazer(DR_effhope_m1,DR_effhope_m3, DR_effhope_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_effhope_m1[,"Std. Error"],HC_DR_effhope_m3[,"Std. Error"],HC_DR_effhope_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af håb på efficacy",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Håb","Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "effhope.tex",
          label="effhope",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")


```

```{r Efficacy contingent on hope}

DR_hopeeff_m1 <- lm(efficacy_index ~ treatment, data = drpanel)
DR_hopeeff_m2 <- lm(efficacy_index ~ treatment + age + male + republican +eduyears + income, data = drpanel)
DR_hopeeff_m3 <- lm(efficacy_index ~ treatment + age + male + republican + eduyears, data = drpanel)

HC_DR_hopeeff_m1 <- coeftest(DR_hopeeff_m1,vcov=vcovHC(DR_hopeeff_m1))
HC_DR_hopeeff_m2 <- coeftest(DR_hopeeff_m2,vcov=vcovHC(DR_hopeeff_m2))
HC_DR_hopeeff_m3 <- coeftest(DR_hopeeff_m3,vcov=vcovHC(DR_hopeeff_m3))


stargazer(DR_hopeeff_m1,DR_hopeeff_m3, DR_hopeeff_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hopeeff_m1[,"Std. Error"],HC_DR_hopeeff_m3[,"Std. Error"],HC_DR_hopeeff_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af håb på efficacy",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid", "Kontrol + Håb","Geografisk + håb","Social + Håb", "Hypotetisk + Håb", "Tid + Håb", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "hopeeff.tex",
          label="hopeeff",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")


```

```{r run interaction efficacy}

stargazer(efficacy_m1, efficacy_m2,efficacy_m3,type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for uddannelse er 'Less than High School", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_efficacy_m1[,"Std. Error"],HC_efficacy_m2[,"Std. Error"],HC_efficacy_m3[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af psykologiske afstande på efficacy betinget af partitilhørsforhold, MTurk",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Hypotetisk", "Social","Tidsmæssig")), c("Alder", "Køn, Mand", "Republikaner")), levels(MTurk$eduOrdered)[-1]),c("Geografisk × Republikaner","Hypotetisk × Republikaner","Social × Republikaner", "Tid × Republikaner")),
          out = "partyeff.tex",
          label="partyeff",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt",
          column.labels=c("(3)", "(4)", "(5)"),
          model.numbers = F)
          
```

```{r Visualise comparison}

#This function is equivalent to dydx in Stata
efficacy_margs <- margins(efficacy_m3, at = list(republican = c(0,1)))
dfmargs1 <- data.frame(summary(efficacy_margs))

ggmargs_polsup_hyp <- ggplot(dfmargs1[dfmargs1$factor %in% "StimulusHypothetical",], aes(x = republican %>% as.factor, y = AME)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), linetype = 3) +
  geom_hline(yintercept=0) + theme_minimal() + scale_x_discrete(labels=c("Ikke-republikanere", "Republikanere")) + labs(x = "", y = "", caption = "95 %'s konfidensintervaller", title = "") + theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank()) + coord_cartesian(ylim = c(-.25,.25))

tikz(file = "gg_polsup_hyp.tex",
     standAlone=F,
     width = 4,
     height = 4,
     sanitize =T,
     )

print(ggmargs_polsup_hyp)

dev.off()

```

```{r Effect of hope and treatment on concern}

DR_hopecon_m1 <- lm(concern_index ~ treatment, data = drpanel)
DR_hopecon_m2 <- lm(concern_index ~ treatment + age + male + republican +eduyears + income, data = drpanel)
DR_hopecon_m3 <- lm(concern_index ~ treatment + age + male + republican + eduyears, data = drpanel)

DR_hopecon_m1div <- lm(concern_index ~ treatment5*treatment_hope, data = drpanel)

HC_DR_hopecon_m1 <- coeftest(DR_hopecon_m1,vcov=vcovHC(DR_hopecon_m1))
HC_DR_hopecon_m2 <- coeftest(DR_hopecon_m2,vcov=vcovHC(DR_hopecon_m2))
HC_DR_hopecon_m3 <- coeftest(DR_hopecon_m3,vcov=vcovHC(DR_hopecon_m3))

stargazer(DR_hopecon_m1,DR_hopecon_m3, DR_hopecon_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hopecon_m1[,"Std. Error"],HC_DR_hopecon_m3[,"Std. Error"],HC_DR_hopecon_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af stimuli på bekymring",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid", "Kontrol + Håb","Geografisk + håb","Social + Håb", "Hypotetisk + Håb", "Tid + Håb", "Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "hopecon.tex",
          label="hopecon",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")


```

```{r Polsup and hope}

DR_hopeints_m1 <- lm(intention_index ~ treatment5, data = drpanel)
DR_hopeints_m2 <- lm(intention_index ~ treatment5*treatment_hope_num + age + male + republican +eduyears + income, data = drpanel)
DR_hopeints_m3 <- lm(intention_index ~ treatment5*treatment_hope_num, data = drpanel)

HC_DR_hopeints_m1 <- coeftest(DR_hopeints_m1,vcov=vcovHC(DR_hopeints_m1))
HC_DR_hopeints_m2 <- coeftest(DR_hopeints_m2,vcov=vcovHC(DR_hopeints_m2))
HC_DR_hopeints_m3 <- coeftest(DR_hopeints_m3,vcov=vcovHC(DR_hopeints_m3))


stargazer(DR_hopeints_m1,DR_hopeints_m3, DR_hopeints_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hopeints_m1[,"Std. Error"],HC_DR_hopeints_m3[,"Std. Error"],HC_DR_hopeints_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af stimuli på handlingsintentioner betinget af håb",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Handlingsintentioner",
          dep.var.labels = "Handlingsintentioner",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid","Ingen stimulus", "Håb")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), c("Geografisk × Håb","Social × Håb", "Hypotetisk × Håb","Tid × Håb")),
          out = "hopeints.tex",
          label="hopeints",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

```

```{r Action and hope}

DR_hopepolsup_m1 <- lm(support_index ~ treatment5, data = drpanel)
DR_hopepolsup_m2 <- lm(support_index ~ treatment5*treatment_hope_num + age + male + republican +eduyears + income, data = drpanel)
DR_hopepolsup_m3 <- lm(support_index ~ treatment5*treatment_hope_num, data = drpanel)

HC_DR_hopepolsup_m1 <- coeftest(DR_hopepolsup_m1,vcov=vcovHC(DR_hopepolsup_m1))
HC_DR_hopepolsup_m2 <- coeftest(DR_hopepolsup_m2,vcov=vcovHC(DR_hopepolsup_m2))
HC_DR_hopepolsup_m3 <- coeftest(DR_hopepolsup_m3,vcov=vcovHC(DR_hopepolsup_m3))


stargazer(DR_hopepolsup_m1,DR_hopepolsup_m3, DR_hopepolsup_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hopepolsup_m1[,"Std. Error"],HC_DR_hopepolsup_m3[,"Std. Error"],HC_DR_hopepolsup_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af stimuli på støtte til politisk handling betinget af håb",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Støtte til politisk handling",
          dep.var.labels = "Støtte til politisk handling",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid","Ingen stimulus", "Håb")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), c("Geografisk × Håb","Social × Håb", "Hypotetisk × Håb","Tid × Håb")),
          out = "hopepolsup.tex",
          label="hopepolsup",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

```

```{r Concern and hope}

DR_hopecon_m1 <- lm(concern_index ~ treatment5, data = drpanel)
DR_hopecon_m2 <- lm(concern_index ~ treatment5*treatment_hope_num + age + male + republican +eduyears + income, data = drpanel)
DR_hopecon_m3 <- lm(concern_index ~ treatment5*treatment_hope_num, data = drpanel)

HC_DR_hopecon_m1 <- coeftest(DR_hopecon_m1,vcov=vcovHC(DR_hopecon_m1))
HC_DR_hopecon_m2 <- coeftest(DR_hopecon_m2,vcov=vcovHC(DR_hopecon_m2))
HC_DR_hopecon_m3 <- coeftest(DR_hopecon_m3,vcov=vcovHC(DR_hopecon_m3))


stargazer(DR_hopecon_m1,DR_hopecon_m3, DR_hopecon_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hopecon_m1[,"Std. Error"],HC_DR_hopecon_m3[,"Std. Error"],HC_DR_hopecon_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af stimuli på bekymring betinget af håb",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid","Ingen stimulus", "Håb")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), c("Geografisk × Håb","Social × Håb", "Hypotetisk × Håb","Tid × Håb")),
          out = "hopecon.tex",
          label="hopecon",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

```

```{r Efficacy and hope}

DR_effhope_m1 <- lm(efficacy_index ~ treatment5, data = drpanel)
DR_effhope_m2 <- lm(efficacy_index ~ treatment5*treatment_hope_num + age + male + republican +eduyears + income, data = drpanel)
DR_effhope_m3 <- lm(efficacy_index ~ treatment5*treatment_hope_num, data = drpanel)

HC_DR_effhope_m1 <- coeftest(DR_effhope_m1,vcov=vcovHC(DR_effhope_m1))
HC_DR_effhope_m2 <- coeftest(DR_effhope_m2,vcov=vcovHC(DR_effhope_m2))
HC_DR_effhope_m3 <- coeftest(DR_effhope_m3,vcov=vcovHC(DR_effhope_m3))


stargazer(DR_effhope_m1,DR_effhope_m3, DR_effhope_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_effhope_m1[,"Std. Error"],HC_DR_effhope_m3[,"Std. Error"],HC_DR_effhope_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af stimuli på efficacy betinget af håb",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Efficacy",
          dep.var.labels = "Efficacy",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Geografisk","Social", "Hypotetisk","Tid","Ingen stimulus", "Håb")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), c("Geografisk × Håb","Social × Håb", "Hypotetisk × Håb","Tid × Håb")),
          out = "effhope.tex",
          label="effhope",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")

```

```{r Hope no comparison}

DR_hope_m1 <- lm(hope_index ~ treatment_hope, data = drpanel)
DR_hope_m2 <- lm(hope_index ~ treatment_hope + age + male + republican +eduyears + income, data = drpanel)
DR_hope_m3 <- lm(hope_index ~ treatment_hope + age + male + republican +eduyears, data = drpanel)

HC_DR_hope_m1 <- coeftest(DR_hope_m1,vcov=vcovHC(DR_hope_m1))
HC_DR_hope_m2 <- coeftest(DR_hope_m2,vcov=vcovHC(DR_hope_m2))
HC_DR_hope_m3 <- coeftest(DR_hope_m3,vcov=vcovHC(DR_hope_m3))


stargazer(DR_hope_m1,DR_hope_m3,DR_hope_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hope_m1[,"Std. Error"],HC_DR_hope_m3[,"Std. Error"],HC_DR_hope_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af intenderet håb på håbsfølelse",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Håb",
          dep.var.labels = "Håb",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Håb","Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "hopehope.tex",
          label="hopehope",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")


```

```{r Hope on concern}

DR_conhope_m1 <- lm(concern_index ~ treatment_hope, data = drpanel)
DR_conhope_m2 <- lm(concern_index ~ treatment_hope + age + male + republican +eduyears + income, data = drpanel)
DR_conhope_m3 <- lm(concern_index ~ treatment_hope + age + male + republican +eduyears, data = drpanel)

HC_DR_conhope_m1 <- coeftest(DR_conhope_m1,vcov=vcovHC(DR_conhope_m1))
HC_DR_conhope_m2 <- coeftest(DR_conhope_m2,vcov=vcovHC(DR_conhope_m2))
HC_DR_conhope_m3 <- coeftest(DR_conhope_m3,vcov=vcovHC(DR_conhope_m3))


stargazer(DR_conhope_m1,DR_conhope_m3,DR_conhope_m2,type="latex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_hope_m1[,"Std. Error"],HC_DR_hope_m3[,"Std. Error"],HC_DR_hope_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af intenderet håb på bekymring",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Bekymring",
          dep.var.labels = "Bekymring",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Håb","Ingen stimulus")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "conhope.tex",
          label="conhope",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r Fear no comparison}

DR_fear_m1 <- lm(fear_index ~ treatment_hope, data = drpanel)
DR_fear_m2 <- lm(fear_index ~ treatment_hope + age + male + republican +eduyears + income, data = drpanel)
DR_fear_m3 <- lm(fear_index ~ treatment_hope + age + male + republican +eduyears, data = drpanel)

HC_DR_fear_m1 <- coeftest(DR_fear_m1,vcov=vcovHC(DR_fear_m1))
HC_DR_fear_m2 <- coeftest(DR_fear_m2,vcov=vcovHC(DR_fear_m2))
HC_DR_fear_m3 <- coeftest(DR_fear_m3,vcov=vcovHC(DR_fear_m3))


stargazer(DR_fear_m1,DR_fear_m3,DR_fear_m2,type="text",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Referencegruppe for indkomst er 'under 100.000'", "Heterogenitetsrobuste standardfejl"),
          se = list(HC_DR_fear_m1[,"Std. Error"],HC_DR_fear_m3[,"Std. Error"],HC_DR_fear_m2[,"Std. Error"]),
          notes.append=FALSE,
          title = "Effekt af intenderet frygt på frygt",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Frygt",
          dep.var.labels = "Frygt",
          dep.var.labels.include = F,
          covariate.labels = append(append(append(append(c("Konstant"),c("Håb")), c("Alder", "Køn, Mand", "Højreorienteret", "Uddannelse, år")),levels(drpanel$income)[-1]), levels(MTurk$eduOrdered)[-1]),
          out = "fearfear.tex",
          label="fearfear",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r Balance tables}

Balance1 <- MTurk %>%
  group_by(Stimulus) %>%
  summarize(AvgAge = mean(age, na.rm = T),
            SDAge = sd(age, na.rm = T),
            minAge = min(age, na.rm = T),
            maxAge = max(age, na.rm = T),
            
            PctWomen = mean(male %>% as.numeric*(-1)+1, na.rm = T),
            
            PctRepublican = mean(republican, na.rm = T),
            sdRepublican = sd(republican, na.rm = T),
            
            ) %>% cbind(MTurk %>% count(Stimulus))

#Count by education
educount <- MTurk %>%
    group_by(eduOrdered) %>% count(Stimulus)

educount_wide <- spread(educount, 
       key = Stimulus,
       value = n)

educount_wide$eduOrdered <-str_replace_na(educount_wide$eduOrdered, replacement = "NA")
educount_wide <- educount_wide %>% remove_rownames %>% column_to_rownames(var="eduOrdered")

#Long to wide
Balance1 <- Balance1[-c(9)]
spread(Balance1, 
       key = Stimulus,
       value = n)

balance_out_wide <- data.frame(Balance1 %>% t)

colnames(balance_out_wide) <- balance_out_wide[1,]
balance_out_wide <- balance_out_wide[-c(1),]

educount_wide <- replace(educount_wide,is.na(educount_wide),0)

#Divider med totaler for procenter
educount_wide_pct <- educount_wide %>% dplyr::select(names(educount_wide)[1])*100/educount_wide %>% dplyr::select(names(educount_wide)[1]) %>% colSums()

for(i in 2:ncol(educount_wide)){
educount_wide_pct <- cbind(educount_wide_pct, educount_wide %>% dplyr::select(names(educount_wide)[i])*100/educount_wide %>% dplyr::select(names(educount_wide)[i]) %>% colSums())
    
}

educount_wide_pct <- rbind(educount_wide_pct, colSums(educount_wide_pct))
educount_wide_pct <- educount_wide_pct %>% dplyr::select(-6)
balance_out_wide <- balance_out_wide %>% dplyr::select(-6)
#End

tutti <- rbind(balance_out_wide,educount_wide_pct)
tutti[] <- lapply(tutti, as.numeric)

tutti %>% mutate(across(where(is.numeric), ~ round(., 2)))


#Lav balancetabel
print(xtable(tutti, caption="Balancetabel", digits=2, type = "latex"), file = "balance.tex")
kable(tutti,format = "latex", digits = 2, caption="Balancetabel")

```

```{r Comparing to ANES}
######ANES######
#Import ANES
anes_timeseries_2020_stata_20220210 <- read_dta("anes_timeseries_2020_stata_20220210.dta")
ANES <- anes_timeseries_2020_stata_20220210 %>%
    mutate_if(
        haven::is.labelled,
        funs(haven::as_factor(.))
    )

ANES$V201510[ANES$V201510 == "-9. Refused"] <- NA
ANES$V201510[ANES$V201510 == "95. Other {SPECIFY}"] <- NA
ANES$V201510[ANES$V201510 == "-8. Don't know"] <- NA

ANES$edu <- ANES$V201510
ANES$edu[ANES$edu == "4. Associate degree in college - occupational/vocational"] <- "5. Associate degree in college - academic"
ANES$edu[ANES$edu == "7. Master's degree (e.g. MA, MS, MEng, MEd, MSW, MBA)"] <-  "6. Bachelor's degree (e.g. BA, AB, BS)"

#An associate degree takes two years to complete and we will have to combine two categories in order to correspond to my data.

ANES_eduount <- ANES %>% count(edu)
ANES_eduount <- ANES_eduount[-c(7), ]

ANES_eduount$procent <- ANES_eduount$n/(ANES_eduount$n %>% sum())*100

######MTurk######
MTurk$edu_compact <- MTurk$eduOrdered
MTurk$edu_compact[MTurk$edu_compact == "Doctorate"] <-  "Professional degree"


MTurk_eduount <- MTurk %>% count(edu_compact)
MTurk_eduount$pct <- MTurk_eduount$n/(MTurk_eduount$n %>% sum())*100
edubind <- MTurk_eduount %>% dplyr::select(-c("n")) %>% cbind(ANES_eduount %>% dplyr::select(-c("n"))) %>% dplyr::select(-c("edu"))
rownames(edubind) <- edubind$edu_compact
edubind <- edubind %>% dplyr::select(-c("edu_compact"))



#age
anes_timeseries_2020_stata_20220210$V201507x[anes_timeseries_2020_stata_20220210$V201507x == -9] <- NA

#Can be recoded as republican: 
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == -1] <- NA
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 10] <- 0
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 11] <- 1
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 12] <- 0
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 20] <- 0
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 21] <- 1
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 22] <- 0
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 30] <- 0
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 31] <- 1
anes_timeseries_2020_stata_20220210$V201075x[anes_timeseries_2020_stata_20220210$V201075x == 32] <- 0

#Sex
anes_timeseries_2020_stata_20220210$V201600[anes_timeseries_2020_stata_20220210$V201600 == -9] <- NA
anes_timeseries_2020_stata_20220210$V201600[anes_timeseries_2020_stata_20220210$V201600 == 2] <- 0

Balance_anes <- anes_timeseries_2020_stata_20220210 %>%
  summarize(AvgAge = mean(V201507x, na.rm = T),
            SDAge = sd(V201507x, na.rm = T),
            minAge = min(V201507x, na.rm = T),
            maxAge = max(V201507x, na.rm = T),
            
            PctWomen = mean(V201600 %>% as.numeric*(-1)+1, na.rm = T),
            sdWomen = sd(V201600 %>% as.numeric*(-1)+1, na.rm = T),

            PctRepublican = mean(V201075x, na.rm = T),
            sdRepublican = sd(V201075x, na.rm = T),
            
            ) %>% cbind(nrow(anes_timeseries_2020_stata_20220210))

Balance_turk_all <- MTurk %>%
  summarize(AvgAge = mean(age, na.rm = T),
            SDAge = sd(age, na.rm = T),
            minAge = min(age, na.rm = T),
            maxAge = max(age, na.rm = T),
            
            PctWomen = mean(male %>% as.numeric*(-1)+1, na.rm = T),
            sdWomen = sd(male %>% as.numeric*(-1)+1, na.rm = T),
            
            PctRepublican = mean(republican, na.rm = T),
            sdRepublican = sd(republican, na.rm = T),
            
            ) %>% cbind(length(MTurk$Stimulus[!is.na(MTurk$Stimulus)]))

balancedf <- Balance_turk_all %>% t %>% cbind(Balance_anes %>% t) 
balancedf <-  as.data.frame(balancedf)
balancedf <-  rbind(edubind,setNames(balancedf, names(edubind)))

kable(balancedf,format = "latex", digits = 2, caption="MTurk/ANES")


```

```{r, balancetests}
# MTurk$Stimulus <- relevel(MTurk$Stimulus, ref = "Control")

test <- multinom(Stimulus ~ age + male + republican, data = MTurk)
relratio <- exp(coef(test))

stargazer(test, 
          type="latex", 
          coef = list(relratio),
          p.auto = F,
          notes.append=F,
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe: Kontrol"),
          title = "Relative risikoratioer for at blive tildelt gruppe",
          header = F,
          single.row = T,
          intercept.bottom = FALSE, 
          dep.var.caption  = "Stimuli",
          dep.var.labels.include = T,
          dep.var.labels = c("Geografisk","Hypotetisk", "Social", "Tid"),
          covariate.labels = c("Konstant", "Alder, år", "Køn, mand", "Republikaner"),
          out = "MTurkmultinom.tex",
          label="MTurkmultinom",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r, Manipulazione}

mc1_MTurk <- lm(ManipulateCloseness ~ Stimulus, data=MTurk)
mc2_MTurk <- lm(ManipulateSocial ~ Stimulus, data=MTurk)
mc3_MTurk <- lm(ManipulateLikelihood ~ Stimulus, data=MTurk)
mc4_MTurk <- lm(ManipulateTimeframe ~ Stimulus, data=MTurk)

HC_mc1_MTurk <- coeftest(mc1_MTurk,vcov=vcovHC(mc1_MTurk))
HC_mc2_MTurk <- coeftest(mc2_MTurk,vcov=vcovHC(mc2_MTurk))
HC_mc3_MTurk <- coeftest(mc3_MTurk,vcov=vcovHC(mc3_MTurk))
HC_mc4_MTurk <- coeftest(mc4_MTurk,vcov=vcovHC(mc4_MTurk))


stargazer(mc1_MTurk,mc3_MTurk,mc2_MTurk,mc4_MTurk, type="latex", out = "manicheck_MTurk.tex",
          star.char = c("\\dag", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001),
          notes = c("† p<0,1; * p<0,05; ** p<0,01; *** p<0,001", "Referencegruppe for afstande er kontrolgruppen", "Heterogenitetsrobuste standardfejl"),
          notes.append=FALSE,
          title = "Manipulationstjek af intenderede effekter, MTurk",
          header = F,
          single.row = T,
          se = list(HC_mc1_MTurk[,"Std. Error"],HC_mc3_MTurk[,"Std. Error"],HC_mc2_MTurk[,"Std. Error"],HC_mc4_MTurk[,"Std. Error"]),
          intercept.bottom = FALSE,
          dep.var.caption  = "Overensstemelse mellem stimuli og oplevet afstand",
          dep.var.labels = c("Geografisk", "Hypotetisk", "Social", "Tid"),
          dep.var.labels.include = T,
          covariate.labels = c("Konstant", "Geografisk", "Hypotetisk", "Social", "Tid", "Ingen treatment"),
          label="manicheck_MTurk",
          decimal.mark = ",",
          font.size = "footnotesize",
          column.sep.width = "-5pt")
```

```{r Runs a factor analysis of concern}
df_factor_concern <- dplyr::select(MTurk,c("concernOverall","concernSocial","concernDistance","concernHypothetical", "concernTime"))

df_factor_concern <- df_factor_concern[complete.cases(df_factor_concern), ]

factanal.concern <- factanal(df_factor_concern, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.concern)

#Scree
fafitfree <- fa(df_factor_concern,nfactors = ncol(df_factor_concern), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_concern<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Bekymring") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of policy support}
df_factor_concern <- dplyr::select(MTurk,c("hightax","airtravel","electriccars"))

df_factor_concern <- df_factor_concern[complete.cases(df_factor_concern), ]

factanal.concern <- factanal(df_factor_concern, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.concern)

#Scree
fafitfree <- fa(df_factor_concern,nfactors = ncol(df_factor_concern), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_polsup<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Policy support") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of intentions}
df_factor_concern <- dplyr::select(MTurk,c("footprint","transport","vacationing"))

df_factor_concern <- df_factor_concern[complete.cases(df_factor_concern), ]

factanal.concern <- factanal(df_factor_concern, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.concern)

#Scree
fafitfree <- fa(df_factor_concern,nfactors = ncol(df_factor_concern), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_intentions<- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "intentions") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r Runs a factor analysis of efficacy}
df_factor_concern <- dplyr::select(MTurk,c("efficacyOthers","efficacySelf","efficacyExternal"))

df_factor_concern <- df_factor_concern[complete.cases(df_factor_concern), ]

factanal.concern <- factanal(df_factor_concern, factors=1, scores = c("regression"), rotation = "varimax", na.action = "na.omit")
print(factanal.concern)

#Scree
fafitfree <- fa(df_factor_concern,nfactors = ncol(df_factor_concern), rotate = "varimax")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
screeplot_efficacy <- ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Faktorer") +
  ylab("Eigenværdi") +
  labs(title = "Efficacy") +
  geom_hline(yintercept=1, linetype=3, size =1) +
  theme_minimal()+ theme(plot.title = element_text(hjust = 0.5))

#screeplot

```

```{r}
tikz(file = "turkscrees.tex",
     standAlone=F,
     width = 6,
     height = 4,
     sanitize =T,
     )

print(plot_grid(screeplot_concern, screeplot_intentions, screeplot_polsup , screeplot_efficacy, 
          ncol =2))

dev.off()
```

